home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / poly / zsolve.c < prev   
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.1 KB  |  87 lines
  1. /* poly/zsolve.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* zsolve.c - finds the complex roots of  = 0 */
  21.  
  22. #include <config.h>
  23. #include <math.h>
  24. #include <stdlib.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_errno.h>
  27. #include <gsl/gsl_complex.h>
  28. #include <gsl/gsl_poly.h>
  29.  
  30. /* C-style matrix elements */
  31. #define MAT(m,i,j,n) ((m)[(i)*(n) + (j)])
  32.  
  33. /* Fortran-style matrix elements */
  34. #define FMAT(m,i,j,n) ((m)[((i)-1)*(n) + ((j)-1)])
  35.  
  36. #include "companion.c"
  37. #include "balance.c"
  38. /* #include "norm.c" */
  39. #include "qr.c"
  40.  
  41. int
  42. gsl_poly_complex_solve (const double *a, size_t n,
  43.             gsl_poly_complex_workspace * w,
  44.             gsl_complex_packed_ptr z)
  45. {
  46.   int status;
  47.   double *m;
  48.  
  49.   if (n == 0)
  50.     {
  51.       GSL_ERROR ("number of terms must be a positive integer", GSL_EINVAL);
  52.     }
  53.  
  54.   if (n == 1)
  55.     {
  56.       GSL_ERROR ("cannot solve for only one term", GSL_EINVAL);
  57.     }
  58.  
  59.   if (a[n - 1] == 0)
  60.     {
  61.       GSL_ERROR ("leading term of polynomial must be non-zero", GSL_EINVAL) ;
  62.     }
  63.  
  64.   if (w->nc != n - 1)
  65.     {
  66.       GSL_ERROR ("size of workspace does not match polynomial", GSL_EINVAL);
  67.     }
  68.   
  69.   m = w->matrix;
  70.  
  71.   set_companion_matrix (a, n - 1, m);
  72.  
  73.   balance_companion_matrix (m, n - 1);
  74.  
  75.   status = qr_companion (m, n - 1, z);
  76.  
  77.   if (status)
  78.     {
  79.       GSL_ERROR("root solving qr method failed to converge", GSL_EFAILED);
  80.     }
  81.  
  82.   return GSL_SUCCESS;
  83. }
  84.  
  85.  
  86.  
  87.